if (compute_equilibrium_effects==T) { 
  mean_y_diff_by_iteration_dist<-data.frame(shocked_person=NA,shocked_degree=NA,  shocked_closeness=NA, shocked_above_med_closeness=NA, month=NA, dist_from_shocked=NA, partial=NA, general=NA)
sum_y_diff_by_iteration<-data.frame(shocked_person=NA, shocked_degree=NA,  shocked_closeness=NA, shocked_above_med_closeness=NA, month=NA, partial=NA, general=NA)

	attach(parameters_new)
  ## shock is increase by one hour
	shock_delta<-1 
	for (month in 1:n_months)
	{
		
		A_graph_this_month<-graph.adjacency(A_array_months[,,month])
		degree_this_month<-degree(A_graph_this_month)/2
		closeness_this_month<-closeness(A_graph_this_month)
		above_median_closeness<-closeness_this_month>=median(closeness_this_month)
		above_median_closeness<-factor(above_median_closeness, labels=c("below med. closeness", "above med. closeness"))
		for (shocked_person in 1:n_students)
		{
		  tmp<-compute_shock_comparative_statics_df(A_array_months[,,month])(solution$mu_y_model, solution$mu_study_model)(solution$s_array_months_model[,month])(shocked_person_index=shocked_person, shock_delta=shock_delta)
		  
		  ## make table summarizing change in y, initial vs final iteration
		  ## denominate change in GPA points
		  sd_y<-1

		  mean_y_diff_by_iteration_dist_tmp<-with(tmp, tapply(y_diff, list(dist_from_shocked, iteration), mean)) / sd_y
		  mean_y_diff_by_iteration_dist_tmp<-data.frame(shocked_person=shocked_person, shocked_degree=degree_this_month[shocked_person], shocked_closeness=closeness_this_month[shocked_person], shocked_above_med_closeness=above_median_closeness[shocked_person], month=month, dist_from_shocked=as.numeric(rownames(mean_y_diff_by_iteration_dist_tmp)), mean_y_diff_by_iteration_dist_tmp)
		  mean_y_diff_by_iteration_dist<-rbind(mean_y_diff_by_iteration_dist, mean_y_diff_by_iteration_dist_tmp)
		  
      ## then computing total gain, exclude the shocked person's achievement gain
      tmptmp<-subset(tmp, dist_from_shocked>=1)
		  sum_y_diff_by_iteration_tmp<-with(tmptmp, tapply(y_diff, list(iteration), sum)) / sd_y
		  sum_y_diff_by_iteration_tmp<-data.frame(shocked_person=shocked_person, shocked_degree=degree_this_month[shocked_person], shocked_closeness=closeness_this_month[shocked_person], shocked_above_med_closeness=above_median_closeness[shocked_person], month=month, t(sum_y_diff_by_iteration_tmp))
		  sum_y_diff_by_iteration<-rbind(sum_y_diff_by_iteration, sum_y_diff_by_iteration_tmp)
		}
	}

	detach(parameters_new)

	rm(mean_y_diff_by_iteration_dist_tmp)
	rm(sum_y_diff_by_iteration_tmp)
	mean_y_diff_by_iteration_dist<-mean_y_diff_by_iteration_dist[-1,]
	sum_y_diff_by_iteration<-sum_y_diff_by_iteration[-1,]
	save(mean_y_diff_by_iteration_dist, file="./output/mean_y_diff_by_iteration_dist")
	save(sum_y_diff_by_iteration, file="./output/sum_y_diff_by_iteration")
} else {
  load("./output/mean_y_diff_by_iteration_dist")
  load("./output/sum_y_diff_by_iteration")
}


melted_mean_y_diff<-melt(mean_y_diff_by_iteration_dist, measure.vars=c("partial", "general"))
names(melted_mean_y_diff)[7]<-"equilibrium"

cutoff_IRF<-4

  ## TABLE 5
  mean_y_diff_table<-data.frame(with(subset(melted_mean_y_diff, dist_from_shocked>=0 & dist_from_shocked<=cutoff_IRF), tapply(value, list(equilibrium, dist_from_shocked), mean)))
  names(mean_y_diff_table)[1:(cutoff_IRF+1)]<-paste(0:cutoff_IRF)
  print(xtable(mean_y_diff_table, label="tab:IRF", caption="Avg. change in GPA by distance from shocked node", digits=3), file="./output/IRF_table.tex")

melted_sum_y_diff<-melt(sum_y_diff_by_iteration, measure.vars=c("partial", "general"))
names(melted_sum_y_diff)[6]<-"equilibrium"

  melted_sum_y_diff<-melted_sum_y_diff[order(melted_sum_y_diff$month,melted_sum_y_diff$shocked_closeness),]
  tmp<-length(order(melted_sum_y_diff$shocked_closeness))/n_months
  tmptmp<-rep((1:tmp)/tmp,2)
  melted_sum_y_diff$shocked_closeness_pctile_by_sem<-tmptmp

  x_lim_tmp<-c(-0.05,1.05)
  y_lim_tmp<-c(0,max(melted_sum_y_diff$value)*1.05)

  ## FIGURE A8, left panel (equilibrium impact of shocks)
  tikz('./output/plot_total_y_diff_by_closeness_pctile.tex', standAlone=F)
    ggplot(data=subset(melted_sum_y_diff, equilibrium=="general"), aes(x=shocked_closeness_pctile_by_sem, y=value, size=shocked_degree, color=equilibrium)) + geom_point(alpha=0.4) + labs(y="total increase in GPA", x="centrality percentile (1 means most central)", size="deg. shocked node") + scale_colour_manual(values = c("blue", "blue")) + scale_shape_manual(values=c(18,16))+my_theme + coord_cartesian(xlim=x_lim_tmp, ylim=y_lim_tmp)
  dev.off()

  ## FIGURE A8, right panel (partial equilibrium impact of shocks)
  tikz('./output/plot_total_y_diff_by_closeness_pctile_partial.tex', standAlone=F)
    ggplot(data=subset(melted_sum_y_diff, equilibrium=="partial"), aes(x=shocked_closeness_pctile_by_sem, y=value, size=shocked_degree, color=equilibrium)) + geom_point(alpha=0.4) + labs(y="total increase in GPA", x="centrality percentile (1 means most central)", size="deg. shocked node") + scale_colour_manual(values = c("red", "red")) + scale_shape_manual(values=c(18,16))+my_theme + coord_cartesian(xlim=x_lim_tmp, ylim=y_lim_tmp)
  dev.off()

source("plot_network_IRF.R")
